Conditions for the cosmological viability of f(R) dark energy models 
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We derive the conditions under which dark energy models whose Lagrangian densities / are written 
in terms of the Ricci scalar R are cosmologically viable. We show that the cosmological behavior 
of f(R) models can be understood by a geometrical approach consisting in studying the m(r) curve 
on the (r, m) plane, where m = Rf,Rn/f t R and r = — Rf,n/f with f t R = Af/AR. This allows 
us to classify the f(R) models into four general classes, depending on the existence of a standard 
matter epoch and on the final accelerated stage. The existence of a viable matter dominated epoch 
prior to a late-time acceleration requires that the variable m satisfies the conditions m(r) ~ +0 
and Am/Ar > -1 at r w -1. For the existence of a viable late-time acceleration we require instead 
either (i) m = — r — 1 , (\^3 — l)/2 < m < 1 and Am/Ar < —1 or (ii) 0<m<latr = —2. 
These conditions identify two regions in the (r, m) space, one for the matter era and the other 
for the acceleration. Only models with a m(r) curve that connects these regions and satisfy the 
requirements above lead to an acceptable cosmology. The models of the type f{R) = aR~ n and 
/ = R + aR~ n do not satisfy these conditions for any n > and n < — 1 and are thus cosmologically 
unacceptable. Similar conclusions can be reached for many other examples discussed in the text. 
In most cases the standard matter era is replaced by a cosmic expansion with scale factor a oc t 1 ^ 2 . 
We also find that f(R) models can have a strongly phantom attractor but in this case there is no 
acceptable matter era. 



I. INTRODUCTION 

The late-time accelerated expansion of the universe is a major challenge to present-day cosmology (see Refs. [3, 0| 
for review). A consistent picture, the concordance model, seems to emerge from the bulk of observations probing the 
background evolution of the universe as well as its inhomogeneities: Supernovae la 0], Cosmic Microwave Background 
anisotropies (CMB) 0], Large Scale Structure formation (LSS) :5|, baryon oscillations weak lensing Q, etc. If 
one assumes today a flat universe with a cosmological constant A and with pressureless matter, observations suggest 
the following cosmological parameters f^A.o ~ 0.7, f2 mi o ~ 0.3 where £lx,o = Px.o/pcr,o for any component X, where 
the subscript stands for present-day values and p cr is the critical density of the universe. 

A cosmological constant term is the simplest possibility to explain the observational data. In fact the recent data 
analysis Q combining the SNLS data with CMB, LSS and the Lyman-cv forest shows that, assuming wbe is constant 
the equation of state parameter of Dark Energy (DE) is found to be wde = — 1.04±0.06 and therefore consistent with 
a cosmological constant. However a cosmological constant suffers from an extreme fine-tuning problem of its energy 
scale if it originates from vacuum energy. For this reason several works have explored alternative explanations, i.e., 
dynamical forms of dark energy. In the absence of any compelling dynamical dark energy model, further insight can 
be gained by considering general models with constant equation of state or some fiducial parametrization, see e.g. 
[Iol |. and [ljj for a recent review. 

The first alternative possibility to a cosmological constant is a minimally coupled scalar field (j>, usually called 
quintessence [H[ • In analogy with inflationary scenarios, this scalar field would be responsible for a stage of accelerated 
expansion while in contrast to inflation this stage occurs in the late-time evolution of the universe. The energy density 
of the scalar field should therefore come to dominate over other components in the universe only recently. This is the 
so called cosmic coincidence problem faced by most dark energy models. In order to alleviate this problem various 
generalizations have been considered, like coupled quintessence models [l3[ in which matter and dark energy scale in 
the same way with time during some epochs. It is however still a challenging task to construct viable scaling models 
which give rise to a matter-dominated era followed by an accelerated scaling attractor flij . 

An important limitation of standard quintessence models is that they do not allow for a phantom regime with 
hide < — 1. A phantom regime is allowed by observations and even favored by some analysis of the data [l5l |. To 
achieve wde < — lj the scalar field should be endowed with a generalized kinetic term, for instance one with a sign 



opposite to the canonical one [1 61 ] . This intriguing possibility is however plagued by quantum instabilities [171 ] . A 
further interesting possibility is provided by non-minimally coupled scalar fields [181 ] and scalar-tensor cosmo logy 
Scalar-tensor DE models can have a consistent phantom regime and a modified growth rate of structure [2l| , 
see also [22[ for a systematic study of the low redshift structure of such theories including a detailed analysis of the 
possibility to have a phantom regime and the constraints from local gravity tests, and [23| for some concrete examples 
of this scenario. In scalar-tensor DE models gravity is modified by an additional dynamical degree of freedom, the 
scalar partner of the graviton. 

Recently there has been a burst of activity dealing with so-called modified gravity DE models (see Ref. [2| for recent 
review and references therein). In these theories one modifies the laws of gravity whereby a late-time accelerated 
expansion is produced without recourse to a DE component, a fact which renders these models very attractive. In 
some models one can have in addition a phantom regime, which might constitute an interesting feature. 

The simplest family of modified gravity DE models is obtained by replacing the Ricci scalar R in the usual Hilbcrt- 
Einstein Lagrangian density for some function f(R). In the first models proposed in DE literature, where a term 1/R 
is added to R [24 , l25l ] , one typically expects that as the universe expands the inverse curvature term will dominate and 
produce the desired late-time accelerated expansion (see Ref. [26[ for a pioneering work in the context of inflation). 
However it was quickly realized that local gravity constraints would make these models non viable [27j (see also 
Ref. (28|). Indeed, it was shown that f(R) models are formally equivalent to scalar-tensor models with a vanishing 
Brans-Dicke parameter wbd = 0. Clearly such models do not pass local gravity (solar system) constraints, in particular 
the post Newtonian parameter 7ppn satisfies 7ppn = 1/2 instead of being very close to 1 as required by observations. 

However, the question of whether local gravity constraints rule out or not f(R) models does not seem to be 
completely settled in the literature [291 ]. Several papers pointed out that local gravity constraints cannot yet rule 
out all possible forms of f(R) theories. For instance, a model containing a particular combination of 1/R and R 2 
terms was suggested (30| and claimed by their authors to pass successfully the solar system constraints, due to a 
large (infinite) effective mass needed to satisfy solar system constraints, and also to produce a late-time accelerated 
expansion (though this latter property does not seem to have been demonstrated in a satisfactory way). Another 
original approach with negative and positive power terms was suggested recently where the positive power term 
would dominate on small scales while the ne gati ve power term dominates on large cosmic scales thereby producing 
the accelerated expansion [3l[ (see however |32|). See Refs. [Hf for a list of recent research in f(R) dark energy 
models. If f(R) models are not ruled out by local gravity constraints it is important to understand their cosmological 
properties. 

Recently three of the present authors [34[ have shown that the large redshift behavior of f(R) = R + aR~ n models 
generically lead to the "wrong" expansion law: indeed, the usual matter era preceding the late-time accelerated stage 
does not have the usual a cx t 2 / 3 behavior but rather a cx t 1 / 2 which would obviously make these models cosmologically 
unacceptable. This intriguing and quite unexpected property of these f(R) models was overlooked in the literature. 
The absence of the standard matter epoch is associated with the fact that in the Einstein frame non-relativistic matter 
is strongly coupled to gravity except for the f(R) theories which have a linear dependence of R (including the ACDM 
model: f (R) = R- A) [p. 

In the Einstein frame the power-law models f(R) cx R~ n (n ^ —1) correspond to a coupled quintessence scenario 
with an exponential potential of a dynamical scalar field. In this case the standard matter era is replaced by a 
"0 matter-dominated epoch" (0MDE) in which the scale factor in the Einstein frame evolves as cx t 3 J 5 [HI]. 
Transforming back to the Jordan frame, this corresponds to a non-standard evolution a oc t 1 / 2 . We wish to stress here 
that cosmological dynamics obtained in the Jordan frame exhibits no difference from the one which is transformed to 
the Einstein frame and transformed back to the original frame. Hence in this paper we shall focus on the analysis in 
the Jordan frame without referring to the Einstein frame. 

This paper is devoted to explaining in detail our previous result and, more importantly, to extend it to all well- 
behaved f(R) Lagrangians. Despite almost thirty years of work on the cosmology of f(R) models, there are in fact 
no general criteria in literature to gauge their validity as alternative cosmological models (see Ref. [35J for one of the 
earliest attempt in this direction). We find the general conditions for a f(R) theory to contain a standard matter era 
followed by an accelerated attractor in a spatially flat, homogeneous and isotropic background. The only conditions 
we assume throughout this paper, beside obviously a well-behaved function f(R) continuous with all its derivatives, 
is that d//di? > 0, to maintain a positive effective gravitational constant in the limit of vanishing higher-order term. 
In some cases however we consider f(R) models which violate this condition in some range of R, but not on the actual 
cosmological trajectories. The main result of this paper is that we are able to show analytically and numerically that 
all f(R) models with an accelerated global attractor belong to one of four classes: 

• Class I : Models of this class possess a peculiar scale factor behavior (a cx t 1 / 2 ) just before the acceleration. 

• Class II : Models of this class have a matter epoch and are asymptotically equivalent to (and hardly distinguish- 
able from) the ACDM model (w c s = — 1). 



• Class III : Models of this class can possess an approximate matter era but this is a transient state which is 
rapidly followed by the final attractor. Technically, the eigenvalues of the matter saddle point diverge and is 
very difficult to find initial conditions that display the approximated matter epoch. 

• Class IV : Models of this class behave in an acceptable way. They possess an approximate standard matter 
epoch followed by a non- phantom acceleration (w c s > — I). 

We can then summarize our findings by saying that f(R) dark energy models are either wrong (Class I), or asymptot- 
ically de-Sitter (Class II), or strongly phantom (Class III) or, finally, standard DE (Class IV). The second and fourth 
classes have some chance to be cosmologically acceptable, but even for these cases it is not an easy task to identify 
the basin of attraction of the acceptable trajectories. We fully specify the conditions under which any given f(R) 
model belongs to one of the classes above and discuss analytically and numerically several examples belonging to all 
classes. 

An important clarification is here in order. It is clear that f(R) gravity models can be perfectly viable in different 
contexts. The most famous example is provided by Starobinsky's model, f(R) = R + aR 2 [3a], which has been the 
first internally consistent inflationary model. In this model, the R 2 term produces an accelerated stage in the early 
universe preceding the usual radiation and matter stages. A late-time acceleration in this model (after the matter 
dominated stage) requires a positive cosmological constant (or some other form of dark energy) in which case the R 2 
term is no longer responsible for the late-time acceleration. 

Our paper is organized in the following way. Section II contains the basic equations in the Jordan frame and 
introduces autonomous equations which are applicable to any forms of f{R). In Sec. Ill we derive fixed points 
together with their stabilities and present the conditions for viable f(R) DE models. In Sec. VI we classify f(R) DE 
models into four classes depending upon the cosmological evolution which gives the late-time acceleration. In Sec. V 
we shall analytically show the cosmological viability for some of the f(R) models by using the conditions found in 
Sec. III. Section VI is devoted to a numerical analysis for a number of f(R) models to confirm the analytical results 
presented in the previous section. Finally we summarize our results in Section VII. We will always work in the Jordan 
frame, in order to show the properties of f(R) models in the most direct way, without the need to convert back from 
the Einstein frame. 



II. f(R) DARK ENERGY MODELS 

A. Definitions and equations 

In this section we derive all basic equations in the Jordan frame (JF), the frame in which observations are performed. 
We will further define all fundamental quantities characterizing our system, in particular the equation of state of our 
system. Actually, as we will see below this is a subtle issue and we have to define what is meant by the Dark Energy 
(DE) equation of state. 

We concentrate on spatially flat Friedman-Lemaitre-Robertson- Walker (FLRW) universes with a time-dependent 
scale factor a(t) and a metric 

d.s 2 = -dt 2 + a 2 (t)dx 2 . (1) 

For this metric the Ricci scalar R is given by 

R = 6(2H 2 + h) , (2) 



(3) 



where H = a/a is the Hubble rate and a dot stands for a derivative with respect to t. 
We start with the following action in the JF 



S = / d 4 xV^ 



1 

2^ 



2 /(i?) + £ rad + £ E 



where k 2 = 8irG while G is a bare gravitational constant, f(R) is some arbitrary function of the Ricci scalar R, and 
£ m and £ ra d are the Lagrangian densities of dust-like matter and radiation respectively. Note that G is typically not 
Newton's gravitational constant measured in the attraction between two test masses in Cavendish-type experiments 
(see e.g. [20|). Then the following equations are obtained (37| 

3FH 2 = K 2 (p m + Prad ) + ^(FR- f)-3HF, (4) 
-2FH = k 2 L m + ^p rad J + F-HF, (5) 



where 

In standard Einstein gravity (/ = R) one has F = 1. In what follows we shall consider the positive-definite forms of 
F to avoid a singularity at F — 0. The densities p m and p ra d satisfy the usual conservation equations 

p m + 3Hp m = , (7) 
p rad + 4Hp Iad = . (8) 

We note that Eqs. ((4]) and (0 are similar to those obtained for scalar-tensor gravity [2l| with a vanishing Brans- 
Dicke parameter uj-qu = and a specific potential U = (FR — /)/2. Note that in scalar-tensor gravity we have 
FR = L so that this term vanishes, while Eq. ^ is similar except for the fact that a kinematic term of the scalar field 
is absent. Hence we can define the DE equation of state in a way similar to that in scalar-tensor theories of gravity 
(see e.g., dHHl]). With a straightforward redefinition of the quantities, we rewrite Eqs. (TJ| and ([5]) as follows 



3F„H 2 = K 2 (pTm + p m + p ra d), (9) 
4 
3' 



-2F H = n 2 ( p m + ~p xad + p DE + pde ) ■ (10) 



We then have the following equalities 

k 2 Pde = \{FR- f)-3HF + 3H 2 (F -F), (11) 
k 2 Pde = F + 2HF-^(FR- f)-(2H + 3H 2 ){F -F). (12) 

The energy density /Ode and the pressure density j>de of DE defined in this way satisfy the usual conservation equation 

Pde = -3H(p U E +Pbe) ■ (13) 

Hence the equation of state parameter wde defined through 

Pde 1 , 2F - 2HF - 4H(F - F) 

wde = = — H i , 14 

Pde (FR- f) - 6HF + 6H 2 (F - F) 

acquires its usual physical meaning, in particular the time evolution of the DE sector is given by 

Pde (z) 



Pde,o 



exp 



, 1 + wde(z') 
Z 1 + z' 



(15) 



where z = ao/a — 1. Note that the subscript "0" stands for present values. It is pde, as defined in Eq. (fTTj) which is 
the quantity extracted from the observations and Wde the corresponding DE equation of state parameter for which 
specific paramctrizations are used. 

Looking at Eqs. ((9|) and (??), one could introduce the cosmological parameters Qx = k 2 Px I (fiF^H 2 ) [H, Ull]. 
However here it turns out to be more convenient to work with the density parameters 

Qx= 3F^' (16) 
where X = m, rad or DE. The quantity wde can further be obtained directly from the observations 

(1 + z) dh 2 /dz - 3h 2 - q. ad ,o(l + zf 
WDE 3[h 2 - a„,o(l + z? - n rad , (l + zY] ' 1 ' 

where h = H/Hq. In the low-redshift region where the contribution of the radiation is negligible, we have 

(l + z)d/i 2 /dz-3/i 2 , , 

wde = p; 7T- — , z < z oq , (18) 

3[h 2 - n m ,o(l + z) A \ 



where z cq is the redshift at which dust and radiation have equal energy densities. Equation (|17p can be extended 
for spatially non-flat universes [39j but we restrict ourselves to spatially flat universes. We also define the effective 
equation of state 



Note that the following equality holds 



2H 

^eff = -1 - ^ ■ (19) 



W e $ = QbE »DE + ^rad , (20) 



if we define & x = n 2 p x /(3F H 2 ). 

B. Autonomous equations 

For a general f(R) model it will be convenient to introduce the following (dimensionless) variables 

*i = ~, (2D 



^ = 7^ = 4 + 2, (23) 



6H 2 H 2 
3FH' 2 



xa = „ TT9 ■ (24) 



From Eq. (01 we have the algebraic identity 



3FH 

It is then straightforward to obtain the following equations of motion 



K 2 

Q m = ~ ~ jj*2 = 1 _ x l — x 2 — %3 — %4 ■ (25) 



where N stands for In a and 



([V = -1 - x 3 ~ 3x 2 + x 1 - xix 3 + x 4 , (26) 

dx2 X\X% , 

dx 4 

-^r = -2X3^4 + Xi X 4 , (29) 



_ dlogF Rf tRR 
m = TT—^ = — i . ( 3 °) 



dlogi? f. R 
dlog/ = Rf,R = X3 
d log R f x 2 



(31) 



where f r = d//di? and J,rr = d 2 f/dR 2 . Deriving R as a function of X3/X2 from Eq. (|3ip . one can express to as 
a function of £3/2:2 and obtain the function m(r). For the power-law model with f(R) = aR~ n the variable to is a 
constant (to = — n — 1) with r = n = 2:3/2:2. In this case the system reduces to a 3-dimensional one with variables 
xi, X2 and 2:4. However for general f(R) gravity models the variable m depends upon r. 
We also make use of these expressions: 

w eS = ~(2x 3 -l), (32) 

1 1-2:4^-2x3 

wde = ~z jz r, (33) 

3 1 - y[l - xi - x 2 - x 3 ) 



where y = F/Fq. 



III. COSMOLOGICAL DYNAMICS OF f(R) GRAVITY MODELS 



In this section we derive the analytical properties of the phase space. 



A. Critical points and stability for a general f(R) 
In the absence of radiation (#4 = 0) the critical points for the system ([26 |) -([28 |) for any m(r) are 



Pi 

P-2 
P 3 
Pi 



(xi,x 2 ,x 3 ) = (0,-1,2), Q m = 0, w eS = -l, (34) 

(ii,a;2,!F3) = (-1,0,0), n m = 2, w eS = l/3, (35) 

(X!,X2,X 3 ) = (1,0,0), tt m = 0, t0eff = l/3, (36) 

(xi, £2,2:3) = (-4,5,0), Q m = 0, w cff = l/3, (37) 



3m 1 + 4m 1 + 4ra \ m(7+10m) m 
P 5 : (xi, £2,0:3)= I— ,-r7^— TT'TTTi— s i "m = 1 -7— ro~, w of f = -— , (38) 



1 + m' 2(1 + m) 2 ' 2(1 + m)/ ' 2(1 + m) 2 ' ° n 1 



m 



/2(l-m) 1 - 4m (1 - 4m)(l + m) \ ^ 2 - 5m - 6m 2 . . 

V 1 + 2m m(l + 2m) m(l + 2m) / 3m(l + 2m) 

where here Q m = 1 — x\ — xi — x-$. 

The points P5 and P§ satisfy the equation x% = — (m(r) + 1)^2, i-e., 

m(r) = -r - 1 . (40) 

When m(r) is not a constant, one must solve this equation. For each root n one gets a point of type P 5 or P 6 with 
m = m(r.i). For instance, the /(-R) = R + aR~ n model corresponds to m(r) = — n(l + r)/r as we will see later, which 
then gives ri 2 = — 1, ti and mi, 2 = 0, —1 — n. If we assume that m = constant then the condition x% = — (m + l)x2 
must hold from Eqs. (|30|) and |3T|) . Hence for m = constant the points i-2,3,5,6 always exist, while P\ and P4 are 
present for m = 1 and m = — 1 respectively. The solutions which give the exact equation of state of a matter era 
(w e g = 0, i.e., a tx i 2/3 or x 3 = 1/2) exist only for m = (P 5 ) or for m = -(5±v / 73)/12 (P 6 ) 0. However the latter 
case corresponds to f2 m = 0, so this does not give a standard matter era dominated by a non-relativistic fluid (4lj |. 

If m (r) is not constant then there can be any number of distinct solutions, although only Pi and those originating 
from P$fi can be accelerated and only P2 and P5 might give rise to matter eras. However P2 corresponds to w e g = 1/3 
and therefore is ruled out as a correct matter era: this is in fact the a cx t 1 / 2 behavior discussed in Ref. [3J| (and 
denoted as 0MDE since it is in fact a field-matter dominated epoch in the Einstein frame). On the contrary, P5 
resembles a standard matter era, but only for m close to 0. Hence a "good" cosmology would be given by any 
trajectory passing near P5 with m close to and landing on an accelerated attractor. Any other behavior would not 
be consistent with observations. 

It is important to realize that the surface X2, £3 for which m(r) = — r — 1 is a subspace of the system (|26M29p and 
therefore it cannot be crossed. This can be seen by using the definition of r and m to derive the following equation 
for r: 

dr (1 , JJ 

=r(l + m + r) , (41) 

AN V J HR' y ' 

which shows explicitly that m = — r — 1 implies dr/dN = as long as R/HR does not diverge. This means that the 
evolution of the system along the m(r) line stops at the roots of the equation m = —r — 1 so that every cosmological 
trajectory is trapped between successive roots. 

In what follows we shall consider the properties of each fixed point in turn. We define m-i = m(P^) and will always 
assume a general m = m(r). 

• (1) Pi: de-Sitter point 
Since w c s = — 1 the point Pi corresponds to de-Sitter solutions (H = 0) and has eigenvalues 



-3, _| ± Y_ZW-, (42) 

where mi = m(r = —2). Hence Pi is stable when < mi < 1 and a saddle point otherwise. Then the condition 
for the stability of the de-Sitter point is given by 

< m(r = -2) < 1 . (43) 



(2) P 2 : </>MDE 

Point P 2 is characterized by a "kinetic" epoch in which matter and field co-exist with constant energy fractions. 
We denote it as a ^-matter dominated epoch (<^>MDE) following Ref. [T^]. The eigenvalues are given by 



m 2 




m 2 



il 



-4 I 12 



m 2 



fr(3 + 4r) 



(44) 



where a prime represents a derivative with respect to r. Hence P 2 is either a saddle or a stable node. If m(r) is 
a constant the eigenvalues reduce to —2,3,4+ l/m2, in which case P 2 is a saddle point. Note that it is stable 
on the subspace X3 = rx 2 for —1/4 < m < 0. However, from Eqs. (f2"T)l and (j2"8)) . one must ensure that the 
term x 3 /m 2 vanishes. Then the necessary and sufficient condition for the existence of the point P 2 is expressed 
simply by 



l im ^1=0. 

212,3^0 m 2 



which amounts to 



f,R 



H 2 f. 



RR 



0, 



(45) 



(46) 



for R/H 2 — ► and f / f^H 2 — > 0. This applies immediately to several models, like e.g., / = log R, R n , R + aR n 
and in general for any well-behaved f(R), i.e., for all the functions that satisfy the condition of application of 
de l'Hopital rule. This shows that the "wrong" matter era is indeed generic to the f(R) models. 

(3) P3: Purely kinetic point 

This also corresponds to a "kinetic" epoch, but it is different from the point P 2 in the sense that the energy 
fraction of the matter vanishes. Point P3 can be regarded as the special case of the point Pq by setting m = 1/4. 
The eigenvalues are given by 



1 



1 



9 + — |r(l + r) ± a i 9 + — |r(l + r) } - 4 <^ 20 + — §r(5 + 4r) 



m 3 



m 3 



m 3 



(47) 



which means that P3 is either a saddle or an unstable node. If m(r) is a constant the eigenvalues reduce to 
2, 5, 4 — l/m3. In this case P3 is unstable for 7713 < and 7713 > 1/4 and a saddle otherwise. 

(4)P 4 

This point has a similar property to P3 because both £l m and w e g are the same as those of P3. It is regarded 
as the special case of the point Pq by setting m = — 1. Point P4 has eigenvalues 



-5, -3, 4(l + l/m 4 ) • 



(48) 



Hence it is stable for —1 < < and a saddle otherwise. Neither of P3 and P4 can be used for the matter- 
dominated epoch nor for the accelerated epoch. 

(5) P5: Scaling solutions 

Point P5 corresponds to scaling solutions which give the constant ratio f2 m /^DE- I n the limit m$ — > 0, it actually 
represents a standard matter era with a oc i 2 / 3 and f2 m = 1. Hence the necessary condition for P5 to exist as 
an exact standard matter era is given by 



m(r = -1) = 0. 



The eigenvalues of P5 are given by 



3(1 + m' 5 ), 



-3rri5 ± \J W5(256r 



160mg 



31m5 — 16) 



4m 5 (m 5 + 1) 

In the limit \m§\ <C 1 the eigenvalues approximately reduce to 

3(l + m' 5 ), ' 



4 V m 5 



(49) 



(50) 



(51) 



The models with 7775 = Rf.Rn/f.R < exhibit the divergence of the eigenvalues as 7775 — > —0, in which case the 
system cannot remain for a long time around the point P5. For example the models f(R) = R — a/ R n with 
n > and a > [24L l25j | fall into this category. An approximate matter era exists also if instead 777.5 is negative 
and non-zero but then the eigenvalues are large and it is difficult to find initial conditions that remain close to 
it for a long time. We shall present such an example in a later section. Therefore generally speaking models 
with 7715 < are not acceptable, except at most for a very narrow range of initial conditions. On the other 
hand, if < 7775 < 0.327 the latter two eigenvalues in Eq. (f50|) are complex with negative real parts. Then, 
provided that 7775 > — 1, the point P5 can be a saddle point with a damped oscillation. Hence in principle the 
universe can evolve toward the point P5 and then leave for the late-time acceleration. Note that the point P2 is 
also generally a saddle point except for some specific cases in which it is stable. Which trajectory (P2 or P5) is 
chosen depends upon initial conditions, so a numerical analysis is necessary. 

Note that from the relation (|4"0")) the condition 777(7-) > is equivalent to r < — 1. Hence the criterion for the 
existence of a saddle matter epoch with a damped oscillation is given by 

777(7- < -1) > 0, m'(r <-!)>-!. (52) 



Note that we also require the condition pi?)) . In order to realize an accelerated stage after the matter era, 
additional conditions are necessary as we will discuss below. Finally, we remark that a special case occurs if 
777 = const. This corresponds to f{R) = —A + aR~ n . In this case the system contains a two-dimensional 
subspace x% = — (m + l)x2 = nx2 and on this subspace the stability of the latter two eigenvalues in Eq. (150[) 
is sufficient to ensure the stability. Working with the {x\, X2, X3) phase space the trajectories that start with 
X3 = nx2, which implies A = 0, remain on the subspace. Then the point is stable in the range < 7775 < 0.327 . 
For A ^ 0, the trajectories start off the subspace and follow the same criteria of stability as for the m ^const. 
case. So there exists a standard saddle matter era for f(R) = —A + aR 1+e with e small and positive. 

• (6) P§: Curvature-dominated point 

This corresponds to the curvature-dominated point whose effective equation of state depends upon the value 
777. It satisfies the condition for acceleration (w c ff < —1/3) when 777,6 < —(1 + \/3)/2, —1/2 < ttiq < and 
777.6 > ( V3 — l)/2. In Fig. [T] we show the behavior of w e s as a function of 777. The eigenvalues are given by 

2 - 37776 - 8ml 2(777|-l)(l+ 777 6 ) 

777 6 ' 7776(1 + 27776) ' 777 6 (1 + 2m 6 ) 

Hence the stability of P6 depends on both 7776 and 777 6 . In the limit 7776 — > ±00 we have P6 — > (—1,0,2) with 
a de-Sitter equation of state (w e g — > —1). This point is stable provided that 777^ > — 1. P6 is also a de-Sitter 
point for 7776 — lj which coincides with Pi and is marginally stable. Since r = —2 in this case, this point is 
characterized by 

m(r = -2)->l. (54) 

It is instructive to see this property in the Einstein frame, i.e. performing a conformal transformation of the 
system [34[. Then one obtains a scalar field with a potential V = (FR — f)/\F\ 2 . This shows that the condition 
TO(r = —2) = 1 corresponds to V.r — 0, i.e., the condition for the existence of a potential minimum. 

The point P6 is both stable and accelerated in four distinct ranges. 
[I] m' & > -1 

When TOg > —1, P6 is stable and accelerated in the following three regions: 

— (A) 7776 < —(1 + v / 3)/2: P6 is accelerated but not a phantom, i.e., w e R > —1. One has w e g — * —1 in the 
limit 7776 ~~ * — 00 • 

— (B) —1/2 < 777 6 < 0: P6 is strongly phantom with w c s < —7.6. 

— (C) 7776 > 1 ; P& is slightly phantom with —1.07 < w e ff < — 1. One has w c g — > — 1 in the limit 7776 — * +°° 
and 7776 - * !• 



[II] 777^ < -1 

When mg < — 1, the point P6 is stable and accelerated in the region 
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Figure 1: The effective equation of state w e g for P§ as a function of m. The point is stable and accelerated in the grayed regions. 
In the region (A) m < — (y/3 + l)/2 the point is always a non-phantom (w c g > —1); in the region (B) —1/2 < m < it is 
strongly a phantom (ui e ff < —7.6); in the region (C) m > 1 it is slightly phantom (—1.07 < w c g < —1) and in the region (D) is 
a non-phantom (ui c ff > —1). In all the other regions Pq is either decelerated or unstable. Notice the gap between w c ff = —1.07 
and —7.6. 

- (D) (V3 - l)/2 < m 6 < 1: here P 6 is a non-phantom, w e s > — 1. 

Therefore from this we derive the first general conclusion concerning f(R) models: the asymptotic acceleration 
cannot have an equation of state in the range —7.6 < w c ff < —1.07. 

If one considers radiation in addition to £1,2.3, then all the points Pi-6 remain the same (with x± = 0) and one 
obtains two additional points 

• (7) P 7 : (x u x 2 ,x 3 ,x 4 ) = (0,0,0,1), O m = 0, w eS = l/3, (55) 

( Am 2m 2m 1 — 2m — 5m 2 \ 1 — 3m 

• 8 P s : (x 1 ,x 2 ,xa,x i ) = I — ,-77- , — 7— — , S2 m = , w eS = — — . (56) 

\1 + m (1 + m) 1 1 + m (1 + m) A J 3 + 3m 

We see that P7 is a standard radiation point. When m(r) is a constant the eigenvalues of P7 are given by 1, 4, 4, — 1, 
which means that P-j is a saddle in this case. The point Ps is a new radiation era (we call it "0 -radiation dominated 
epoch" ) which contains non-zero dark energy. Since the effective equation of state is constrained by nucleosynthesis 
to be close to 1/3, Ps is acceptable as a radiation epoch only for mg close to 0. 

The eigenvalues of Ps are given by 

1 4 /, +m M m 8 - 1 ± y/81ml + 30m 8 - 15" 

1, 4(1 + m 8 ), 2(71*8 + 1) • (57) 

In the limit mg — > the last two are complex with negative real parts, which then shows that Pg is a saddle around 
the radiation point. Hence the solutions eventually repel away from the radiation era and are followed by one of the 
fixed points given above. Unlike the matter point P5 there are no singularities for the eigenvalues of Pg in the limit 
mg — * 0. We also note that Pg is on the line m = — r — 1 as in the case of the matter point P5. If the condition for the 
existence of the matter point P5 is satisfied (i.e., m ~ and r ~ —1), there exists a radiation point Pg in the same 
region. Then a viable cosmological trajectory starts around the radiation point Pg with m w and then connects to 
the matter point P5 with m « 0. Finally the solutions approach either of the accelerated points mentioned above. 

IV. FOUR CLASSES OF MODELS 




For a cosmological model to work, it has to possess a matter dominated epoch followed by an accelerated expansion. 
In our scenarios this would be a stable acceleration (late-time attractor). We require that the matter era is long enough 
to allow for structure formation and that an effective equation of state is close to w e g = in order to match the 



observations of the diameter distance of acoustic peaks of CMB anisotropies, i.e., it has to expand as a ~ i 2 / 3 . Now 
we study the conditions under which these requirements are met. 

Let us recall again that P2 exists as a saddle or a stable node. Then the 0MDE is always present provided that 
the condition (|45|) is satisfied and only by a choice of initial conditions one can escape it. Hence below we examine 
the cases in which initial conditions exist such that the standard matter era P5 for \m\ 1 is also a saddle. When 
this is possible, a numerical analysis is necessary to ascertain the basin of attraction of P2 and P5. In particular, it is 
necessary to see whether initial conditions that allow for a radiation epoch lead to Pi or P5. Then, if P5 exists and 
is a saddle, we examine the conditions for a late-time accelerated attractor. 

A. Transition from the matter point P5 to an accelerated point P§ or Pi 

The only point which allows for a standard matter era is P5 when m(— 1) — > 0, so this is the first condition for 
a theory to be acceptable. If m(— 1) is non-vanishing the matter epoch can be characterized by a ~ ^2(i+m)/3^ 
which is still acceptable if |m| <C 1. So from now on when we write m(— 1) — > we always mean |m(— 1)| <C 1. 
The corresponding point P5 with |mg| <C 1 will be denoted as P 5 . In the general case, Eq. (|40j) has several roots 
^a.fc.. and therefore m a j,_ and correspondingly there will be several points Pg( a ,i>,...), ^6(0,6,..) ■ Let us call the line 
m = —r — 1 on the (r, m) plane the critical line, since the points P5 and Pg lie on this line. From the matter epoch 

Pg at (r,m) — (—1,0), the trajectories can reach an acceleration point at either Pi or one of the points P5 (beside 

P^) or Pg, the only points that can be accelerated. The point Pi is stable and accelerated only for < mi < 1. The 
point P5 corresponds to an accelerated solution for 777,5 > 1/2 and 7775 < — 1; however, it can be shown that it is not 
stable (saddle or unstable node) in both regions. Therefore wc only need to study the transition from the matter point 

P5 to the accelerated point Pg. Generally speaking, a f(R) model is cosmologically viable if one of the transition 

P 5 (0) -> P 6 or P 5 (0) -> Pi is possible. 

The point Pg : (r, 777) = (—1,0) can be approached from the positive 777 side or from the negative one. In the 
first case, two eigenvalues are complex while the real part of the eigenvalues remains finite and negative. Then the 
trajectory exhibits a damped oscillation around the matter point, before leaving for the acceleration. In the second 
case, the eigenvalues are real and diverge for 777 — > —0. Then the matter era is very short and it is very difficult to 
find initial conditions that lead to a successful cosmology. The pure power-law model f(R) = aR~ n is a special case 
because then Pg is actually stable for m = —1 — 77 small and positive, so it is not possible to reach the acceleration 
at Pg [note that in this case the system is two-dimensional with the latter two eigenvalues in Eq. I|5ip], For the model 
f(R) = —A + aR 1+e with e small and positive, the transition from P^ to Pg is instead possible and these models 
are cosmologically acceptable. This shows that a ACDM cosmology is recovered for this model in the limit e — » +0 
but not e — > —0. 

As we have seen in the previous section, the point Pg is stable and accelerated in four distinct regions: (A) 
777g < -(1 + V3)/2, (B) -1/2 < m 6 < 0, (C) TOg > 1 (all these are stable if m 6 > -1); and finally, if m' 6 < -1, (D) 
— l)/2 < 777g < 1. In the regions (A) and (D) the point Pg leads to a non-phantom acceleration with w c s > — 1, 
whereas the region (B) corresponds to a strongly phantom (w c h < —7.6) and the region (C) to a slightly phantom 
(—1.07 < TTJoff < —1). In what follows we shall discuss each case separately. 

1. From P 5 (7715 > -l,77i > 0) to Pi or to P 6 (mj > -1) in the regions (A), (B), (C) 

In the positive m region the matter point Pg is a saddle for m' 5 > — 1. We require the condition m' 6 > —1 for the 
stability of the point Pg in the regions (A), (B) and (C). Let us then assume that beside the root at m ~ +0 there 
are three roots which exist in the regions (A), (B), (C), i.e. mg a , mgb, 777g c , respectively. A good cosmology goes 

from a saddle Pg to a stable acceleration, either Pg a ,P6&,P6c or Pi. Now P^ is a saddle if m' 5 > — 1, while Pg is 
stable if m' 6 > —1. This shows that the curve m(r) must intersect the critical point line m = —r — 1 with a derivative 
777 5 6 > — 1. If the intersection occurs with a derivative m' 5 6 < —1, the cosmological model is unacceptable, either 
because the matter era is stable or because the accelerated epoch is not stable. 

We can therefore draw on the (r, m) plane the "forbidden direction regions" around the critical points, i.e. the 
direction for a curve m(r) intersecting the line 777 = — r — 1 that must not be realized (see Fig. [3 where we plot several 
possible m(r) that belongs to four general classes as detailed below). So, for any given m(r) model, one has simply 
to look at the intersections of m(r) with m = — r — 1 to decide if that model passes the conditions for a standard 
matter-acceleration sequence. Generally speaking, if the m(r) line connects the standard matter era (r, m) = (—1,0) 
with an accelerated point Pg or Pi without entering the forbidden direction region, then that model is cosmologically 



viable. Otherwise, either because there is no connection at all or because the connection has the wrong direction, 
then the model is to be rejected. 

In general, of course, any m(r) line is possible. However, assuming F > 0, one sees that r(R) is a monotonic 
function, and therefore m(r) is single-valued and non-singular (remember we are assuming a regular / with all its 
derivatives). This simple property is what we need to demonstrate our claims. In fact, it is then simple to realize 
by an inspection of Fig. [2] that indeed it is impossible to connect points near m = +0 with points in (A), (B) or 
(C). To do so it would require in fact either entering the forbidden direction regions or a turn-around of m(r), i.e. a 
multi-valued function, or a singularity of to at finite r, or finally a crossing of the critical line. This simple argument 
shows that the matter era with m w +0 cannot connect to Pq in the region (A), (B) or (C). Hence the only accelerated 
point left is Pi (which is stable only for < m(r — —2) < 1). Notice that this argument applies for any number of 
roots in (A), (B) or (C). 

A connection to Pq is however possible at r — ► ±oo, with slope m' 6 = —1, i.e. when the curve m(r) is asymptotically 
convergent on the m = — r — 1 line. Even in this case, the final acceleration is de-Sitter, although with X2, X3) — 
(—1,0,2) instead of Pi : (xi, X2, X3) = (0,-1,2). To complete this demonstration we need also to ensure that although 
the m(r) line can have any number of intersection with the critical line, no cosmological trajectory can actually cross 
it. This property is indeed guaranteed by Eq. (|41[) : trajectories stop at the intersections of m(r) with the critical line 
and remain trapped between successive roots. 

2. From P5 (m < 0) to Pe {m' 6 > —1) in the region (B) 

There is then a further option: P^ in the (B) region, i.e. 7715 < 0. When m is close to —0, one of the last two 

eigenvalues in Eq. (f5Tj) is positive whereas another is negative. This shows that in this case the point P^ is a saddle 
independently of m' 5 . Note that the accelerated point Pq in the region (B) is stable for m 6 > — 1. 

Let us first consider the case m' 5 > — 1. Then the same argument applies for the positive to case discussed above. 
The m(r) curves can not satisfy both the conditions m' 5 > — 1 and to 6 > — 1 required for the existence of the stable 
accelerated point Pq in the regions (A), (B) and (C). However there is one exception. If the matter root m is small and 
strictly negative and m' 5 > —1, then Pq for the same root lies in the (B) region and is a valid acceleration point. In 
other words, P5 and Pq coincide in the (r, to) plane and are both acceptable since m' 5 6 > — 1. The simplest possibility 
is m = const 6 (—1/2,0). For instance, for the power-law models f(R) — aR 9 (i.e. to = —0.1), the transition 

from an approximately matter epoch P^ to an accelerated era Pq is possible. However the matter period is short 
because of real eigenvalues which diverge in the limit to — —0, Another possibility is to = a + br, i.e. a straight line 
intersecting the critical line at some point with abscissa (a — b)/(l + b) G (—1/2, 0) and a slope b > — 1. 

When to' 5 < —1 it is possible to reach the stable accelerated point Pq in either of the regions (A), (B), (C) with 
to 6 > — 1. However the matter epoch does not last long in this case as well because we have seen an eigenvalue is 
very large. Moreover, by construction there will always be the final attractor Pq in the region (B) for the same to, 
whose effective equation of state corresponds to a strongly phantom (w e g < —7.6). 

Thus if the matter point Pg exists in the region m < 0, the models are hardly compatible with observations because 
the matter era is practically absent and because most trajectories will fall in a unacceptable strongly phantom era. 

3. From P5 (7715 > — l,m > 0) to Pe (m' 6 < —1) in the region (D) 

We come to the fourth range, i.e. the region (D). Now the situation is different for the point Pq, since to 6 has to be 

less than —1 in order to be stable. Then it is possible to leave the matter epoch P^ (which satisfies m' 5 > — 1, to > 0) 
and to enter the accelerated epoch Pq (to 6 < — 1) as we illustrate in Fig. [2] (Class IV panel). Therefore these models are 
compatible with standard cosmology: they have a matter era followed by a non-phantom acceleration with w c g > — 1. 
Note that the saddle matter epoch needs to be sufficiently long for structure formation to occur. Later we shall 
provide an example of such models. 

Finally, we must mention an exception to this general argument. If the m(r) line has a derivative exactly to' = — 1 
at the critical point, then that point is marginally stable and our linearized analysis breaks down. In this case, 
one has to go to a second-order analysis or to a numerical study. We will encounter such a situation for the model 
f(R) — R\og(aR) q we study later. The same applies if to' — * ±00, i.e. for trajectories that lie on the borders of the 
forbidden regions. 




Figure 2: The (r, m) plane for the four classes of f(R) models. In all panels, the straight diagonal line is the critical line 
m = — r — 1. In the dotted ranges Pg is not accelerated or is unstable, if we assume m' 6 > — 1. In the thick ranges labelled by 
A,B and C, P§ is accelerated and stable, again assuming m' 6 > — 1 (we omit the region D for clarity except for the Class IV 
panel). The gray triangles represent the forbidden directions near the critical points. The dashed green lines are hypothetical 
m(r) curves, intersecting the critical line in the critical points P5 and P§. The intersection at (r, m) = (—1,0) (light gray 
triangles) corresponds to the standard matter epoch P^\ In Class I models, the m(r) curve does not intersect (r, m) = ( — 1, 0) 
and therefore there is no standard matter era. In Class II models, the point (r, m) = (—1,0) is connected to the Pi de-Sitter 
point (along the segment 0<m<latr = —2) and therefore represents a viable cosmological solution. The two additional 
critical points in the regions A and C are unstable since the curve enters the forbidden triangles and are therefore not acceptable 
as final accelerated stages. In Class III models the m(r) line with a slope m' > — 1 intersects the critical line at a negative m 
in the strongly phantom range (B). Note that the curves with m' 5 < — 1 which are attracted by P§ in the region (A), (B), (C) 
are possible, but such cases are not viable because of the absence of a prolonged matter era for m < 0. In Class IV models, 
the m(r) curve connects the matter era with m' 5 > — 1 to the region (D) with a derivative m' 6 < — 1 and therefore represents 
a viable cosmology with a matter era followed by a stable acceleration (w e g > — 1). No single trajectory can cross the critical 
line m = — r — 1: each solution is trapped between two successive roots on the critical line. 



B. Classification of f(R) models 

These discussions show that we can classify the f(R) models into four classes, as anticipated in Introduction. The 
classification can be based entirely upon the geometrical properties of the m(r) curve and applies to all the cases in 
which an accelerated attractor exists (see Fig. [SJ. 

• Class I : This class of models covers all cases for which the curve m(r) does not connect the accelerated attractor 
with the standard matter point (r, m) = (—1,0), either because m(r) does not pass near the matter point, i.e. 
m(r — ► —1) 7^ , or because the branch of m(r) that accelerates is not connected to (r,m) = (—1,0). Instead 
of having a standard matter phase, the solutions reach the 0MDE fixed point P2 with a "wrong" evolution of 
the scale factor (a cx t 1 ^ 2 ) or bypass it altogether by falling on the final attractor without a matter epoch at all. 
The final accelerated fixed points, if they exist, can be in any of the three ranges of Pq. 



• Class II : For these models the m(r) curve connects the upper vicinity of the point (r, m) = (—1,0) (with m > 



and m' 5 > — 1) to the point Pi located on the segment < m < 1 at r = —2, or asymptotically to P$(r — > ±oo). 
Since the approach to P5 is on the positive side of m, the trajectory exhibits a damped oscillation around the 
matter point [see Eq. (|51|)]. which is followed by the de-Sitter point Pi or Pq(i- — > ±00). Models of the Class II 
are observationally acceptable and the final acceleration corresponds to a de-Sitter expansion. 

• Class III: For these models the m(r) curve intersects the critical line at— l/2<m<0 (i.e. region B). In all 
these cases the approximated matter era is a very fast transient and only a narrow range of initial conditions may 
allow it. Generically, the matter era is followed by a strongly phantom acceleration, although one could design 
models with the other ranges of the critical line. The closer to a standard matter epoch, the more phantom 
the final acceleration is (u> c ff — * — 00 as m — ► —0). Since the matter era is practically unstable and the highest 
effective equation of state is —w c s = 7.6 (which implies wde — w c ff,o/^DE,o even smaller), these models are 
generally ruled out by observations (although a more careful numerical analysis is required). 

• Class IV: For these models the m(r) curve connects the upper vicinity of the point (r, m) = (—1,0) (with 
m' 5 > —1, m > 0) to the region (D) located on the critical line m = —r — 1 (with m' 6 < —1). These models are 
observationally acceptable and the final acceleration corresponds to a non-phantom effective equation of state 

(Wolf > -1). 

In Fig. [3] we show a gallery of m{r) curves for various f(R) models. The above discussions clarify the conditions for 
which f(R) dark energy models are acceptable. Only the Class II or Class IV models are in principle cosmologically 
viable. However we need to keep in mind that what we have discussed so far corresponds to the behavior only around 
critical points. One cannot exclude the possibility that single trajectories with some special initial conditions happen 
to reproduce an acceptable cosmology. It is therefore necessary to confirm our general analysis with a thorough 
numerical check; by its nature, this check can only be done on a case-by-case basis, and to this we turn our attention 
in the next sections. 



V. SPECIFIC MODELS: ANALYTICAL RESULTS 

In this section we shall consider a number of f(R) models in which m can be explicitly written in terms of the 
function of r and study the possibility to realize the matter era followed by a late-time acceleration. Most of the 
relevant properties of these models can be understood by looking at the m(r) curves of Fig. [31 

1. f(R) = aR~ n 

This power-law model gives a constant m from Eq. (|30p , namely 

m = -n - 1 , (58) 

where r — n. The curve m(r) degenerates therefore to a single point and this case reduces to a two-dimensional 
system in the absence of radiation because of the relation 2; 3 = 11x2- Hence the condition m(r = — 1) = is satisfied 
only for n = — 1, i.e. Einstein gravity. Since the initial conditions around the end of the radiation era are given for 
positive R, the positivity of the term f t R = ~naR~ n ~ 1 requires that a < for n > and a > for n < 0. From 
Eq. ([55)1 one has m > for n < — 1 and m < for n > — 1. Since the eigenvalues of P5 for this model are given by the 
latter two in Eq. ([5U|) [4 if, the matter point P5 is a stable spiral for n < — 1 (around m — » +0). Then the solutions 
do not leave the matter era for the late-time acceleration. 

On the other hand, P5 is a saddle point for — 1 < n < —0.713 while the </>MDE point P2 is stable in the overlapping 
range — 1 < n < —3/4. However one of the eigenvalues of P5 exhibits a positive divergence in the limit m — > — 0, 
which means that the matter point becomes repulsive if m is very close to —0. As we anticipated, in the region 
around m = —0 the effective equation of state for Pq corresponds to the strongly phantom type (w e s < —7.6), i.e., 
to our Class III models. [4l|. The above discussion shows that the saddle point P5 is connected to either the 0MDE 
point P2 or the strongly phantom point Pq. The more one tries to get a standard matter era for n — * — 1, the more 
phantom becomes the final acceleration and the more divergent becomes the eigenvalues. Moreover if we take into 
account radiation, the solutions tend to stay away from the point P5. 

So the models of this type are always in Class I except for (i) — 1 < n < —0.713 (Class III) and for (ii) —1.327 < 
n < — 1 (they are asymptotically not accelerated). Similar conclusions were found in Ref. [43| . 

The pure power-law models correspond to points (r = —n, m = — 1 — n) in the (r, m) plane. We can notice that the 
ACDM model / = R — A corresponds to the horizontal line m = 0, which connects the matter era at (r, m) — (—1,0) 
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Figure 3: This figure illustrates several possible m(r) curves (thick dashed green line). Only the / = R\ogR and the / = 
Rexp(l/R) models show an acceptable connection between the matter point (r, m) = ( — 1,0) and the de-Sitter point Pi along 
the dashed segment at r = —2. In all other cases, there is either no intersection of the m(r) curves with the critical line 
m — —r — 1 near (r, m) = (— 1, 0) or the m(r) curve enters the forbidden direction regions (the gray triangles). In all panels we 
show the forbidden regions for three points in the (A,B,C) ranges of P§, even when there are no critical points in one of those 
regions. For clarity we omit the range (D). 



with the de-Sitter acceleration P\ at (r, m) = (—2, 0) and is therefore a valid Class II model. A possible generalization 
of ACDM is given by the models 

/(i?) = (i? h -A) c , (59) 

which generate a tilted straight line m(r) = r(l — c)/c+b— 1. If the intersection m = —1 + be with the critical line is 
at < m -C 1 and the slope is given by —1 < (1 — c)/c < 0, then the matter era is connected with Pi and the model 
is acceptable (Class II). 



2. f(R) = R + aR~ 



This model was proposed in Refs. [24],[25| to give rise to a late-time acceleration. From Eqs. (f3T)|) and (|31[) we obtain 

m(r) = - n(1 + r) . (60) 
r 

Notice that m(r) is independent of a. Since m(r = — 1) = the models satisfy the necessary condition for the 
existence of the matter point P5 . 

Let us analytically study the attractor behavior of the model in more details. Substituting Eq. (|4T)|) for Eq. (|60[) . 

we find the solutions m a = or = — (n + 1), which holds for the points P5 and Pq. In this case the points P5 and 
Pq are characterized by 

p 5a - (°>-^4)' ° m=i ' ^ff=°' ( 6i ) 



, 3(n + l) 4n + 3 4n + 3\ 8n 2 + 13n + 3 „ 1 

: , 2 ,— ^ , = tt- 2 , W eS = -l--, (62) 

n 2n z 2n I 2n z n 



. 2(n + 2) An + 5 n(An + 5) \ „ 6n 2 + 7n - 1 ,„„, 

p ®> '■ — n — rr>7 — 1 i\/o — r~rT'7 — ttttt; — r~rr > "m = 0, w cff = -— — — — — — . (63) 



2n+l ' (n+ l)(2n + l)' (n+l)(2n+l) J ' ' 3(n + l)(2n+ 1) 

Note that for m a = 0, Pe a goes to infinity. We are interested in the case where a (quasi) matter era is realized around 

7715 w 0. 

This family of models splits into three cases: 1) n < — 1, 2) — 1 < n < 0, and 3) n > 0. The intermediate cases 
n = 0, 1 are of course trivial. 

• Case 1 (n < —1). 

Since m! = n/r 2 we see that m'{— 1) < —1 and therefore the matter epoch around m « +0 is stable and no 
acceleration is found asymptotically (Pi is stable as well for —2 < n < 0). The case n = — 2 corresponds to 
Starobinsky's inflation model and the accelerated phase exists in the asymptotic past rather than in the future. 
This case does not belong to one of our main classes since there is no future acceleration. 

• Case 2 (-K n < 0). 

Then the condition at r = — 1 is fulfilled for R — » 00, and we see that m = n(n + l)aP _,l_1 /(l — naR^ 71 ^ 1 ) 
approaches zero from the positive side if a < 0. In this case, there are damped oscillations around the standard 
matter era and the final stable de-Sitter point Pi can be reached (P$ is unstable): this is the Class II model. 
Notice that F < for small R, but F > along the cosmologically acceptable trajectory. When a > 0, two of 
the eigenvalues diverge as m — > — and the matter era becomes unstable. 

• Case 3 (n > 0). 

In this case the stable accelerated point Pq exists in the non-phantom region (A) because of the condition 
m = —n — 1 < — 1. If a > 0, m approaches zero from the positive side. Then there are oscillations around the 
matter era but the accelerated point Pi is unstable (since mi = —n/2 < 0). Since m' 5 — n > 0, the matter era 
corresponds to a saddle. However P5 with m' 5 > — 1 cannot be connected to Pq in the region (A), as we showed 
in the previous section. Hence we do not have a stable accelerated attractor after the matter epoch. When 
a < 0, m approaches zero on the negative side and here again the matter point becomes effectively unstable 
since one of the eigenvalues exhibits a positive divergence. Then this case does not possess a prolonged matter 
epoch and belongs to the Class I. The first panel of Fig. [3] shows graphically why models like f(R) = R + a/R 
cannot work as a viable cosmological model: the accelerated point is disconnected from the matter point. 

In the next section we shall numerically confirm that the matter phase is in fact absent prior to the accelerated 
expansion except for models f(R) = R + aR~ n with a < and —1 < n < 0. In any case, all these power-law cases are 
cosmologically unacceptable. These results fully confirm the conclusions of Ref. [34| reached by studying the Einstein 
frame. The single exception pointed out above for —1 < n < was not a part of the cases considered in Ref. [34| . 
since F < for small R. 



3. f(R) = R p exp(qR) 



In this model m is given by 

m(r) = — r + — . (64) 
r 

Notice that for the pure exponential case (p = 0) we have to = —r and £3/771 — > x% — > so that P2 exists while P5 
does not. Otherwise the function to vanishes for r — > ±^/p, which means that the condition (|49|) for the existence of 
the matter era holds only for p = 1. However, since in this case m'(r = —1) = — 2 < —1 , the point P5 is a stable 
spiral for m > 0. So the entire family of models is in fact ruled out. 

In the limit m — > +0, Pg can not be used for the late-time acceleration in addition to the fact that P5 is stable. 
Moreover since m(r = —2) = 3/2 for p = 1, the de-Sitter point Pi is not stable. We note that Eqs. (f64|) and (f40|) 
are satisfied in the limit m — > +00 and r — > —00, see Fig. [3J Since the eigenvalues in Eq. (|53p are —4, —4, in this 
case, the point Pe:(xi,X2,Xa) = (—1,0,2) with to — > +00 is marginally stable with an effective equation of state 
w e g — > — 1. In fact, when to > 0, we have numerically checked that the final attractor is either the matter point P5 
or Pq : (x±, X2, Xa) — (— 1, 0, 2) (but then without a preceding matter phase), depending upon initial conditions. 

Thus models of this type do not have the sequence of matter and acceleration for p = 1, whereas the models with 
p / 1 belong to Class I. 



I f (R) = BP (log oRY 

In this model we obtain the relation 

p 2 + 2pr - r(q - r + qr) 

m(r) = . (65) 

qr 

Since m(r = —1) = — (p— l) 2 /g, the matter epoch exists only for p = 1. When p = 1 one has m(r = —2) = 1 — 1/ (2q), 
which means that P± is stable for q > whereas it is not for q < 0. The derivative term m'(r) is given by 

m'(r) = -l + — (66) 

Since m'(r = —1) = —1 the point P5 is marginally stable. However we have to caution that m does not exactly 
become zero. In fact when r < —1 we have m'(r) > —1 and m(r) > for q > 0, which means that the quasi matter 
era with positive to is a saddle point. Similarly the accelerated point Pq in the region (C) is stable for q > whereas 
it is not for q < 0. Hence both Pi and Pq are stable for positive q. However one can show that the function m(r) 
given in Eq. (|65p satisfies m(r) < — r — 1 in the region r < —1 for p = 1 and g > 0. Hence the curve fl55J) does not 
cross the point Pg in the region (C). Then the only possibility is the case in which the trajectories move from the 
quasi matter era P5 to the de-Sitter point Pi. In the next section we shall numerically show that the sequence from 
P5 to Pi is in fact realized. 

Thus when p = 1 and q > the above model corresponds to the Class II, whereas the models with p / 1 are 
categorized as the Class I. 



,5. f(R) = R p exp(q/R) 

This model gives the relation 

p + r(2 + r) 

TO(r) = , (67) 

r 

which is independent of q. Here we have m(r = —1) = p — 1, so a matter era exists for p = 1. In this case one has 
m(r) — — (7- + l) 2 jr > for r < 0. Since m(r = —2) = 1/2 for p = 1, the point Pi is a stable spiral. The derivative 
term m'(r) is given by m'(r) = —1 + 1/r 2 , which then implies m'(r = —1) = and m'(r < —1) > —1. This shows 
that P5 is a saddle whereas Pq in the region (C) is stable. The curve ([FT]) satisfies the relation m(r) < — r — 1 in the 
region r < —1 for p = 1 and also has an asymptotic behavior m(r) — » — r in the limit r — > —00. Then in principle it 
is possible to have the sequence P5 — > Pe(r — ► —00), but the trajectory from the point P5 is trapped by the stable 
de-Sitter point Pi which exists at (r, to) = (—2, 1/2). We note that one of the eigenvalues for the point P5 is large 
(3(1 + m' 5 ) = 3) compared to the model f(R) = R(\ogaR) q whose eigenvalue is close to (but positive) around 
m = 0. In such a case the system does not stay around the matter point P5 for a long time as we will see later. 
Thus the model with p = 1 belongs to the Class II, whereas the models with 1 correspond to the Class I. 



6. f(R) = R + aR 2 - A 



In this case the function m(r) is given by 

m.(r\ = — 

1 + A(r) 

where 



m(r) = ' 1 " r +^ (r) , (68) 



A(r) = + r) 2 + Aar(2 + r) , a = aA. (69) 
Here we assume that a, A > 0. The equation, m(r) = — 1 — r, gives three solutions 

1 + 4a ± 2P 

^ 2 = - l + 4« ' 2 ' (70) 



where B = y/a(l + 45). Then we obtain three points P5 and three Ps. For P5 we have 

/ 6a B(B±8a) 8a±B \ 
n a< b-{x 1 ,x 2 ,Xs)-y 2&±B , 2 (P ± 2a) 2 ' 4a ± 2Py< ' { ' 

P 5 , c :(xi,X2,x 3 )= ■ (72) 

The point Ps )C is unphysical since ri m < 0. The points P$ a ,b reduce to a matter point in the limit d <C 1. At the 
lowest order one has w e ff ~ T4Vfi/3. This shows that a standard matter era can exist either for a — > 0, i.e., for 
the ACDM model, or for A — > 0, i.e., for the Starobinsky's model f(R) = R + aR 2 . In the limit a — > the only 
accelerated point is the de-Sitter point P\. Since the condition m{— 2) = 1 is satisfied for any a, we see that this /(P) 
model is always attracted by the de-Sitter acceleration. 
Models of this type belong to Class II. 



7. f{R) = R- ^\/R + ^ 2 R 2 

This model was proposed in Ref. [30]. In this case Eq. (|40|) reads: 

P 3 + P 2 (1 + r) + $ (1 - r) = , (73) 

where P needs to be real solutions. Since the solutions for this equation are quite complicated, we will not write them 
down here. The necessary condition for the existence of the matter phase is as usual m(— 1) = 0. We have here 

Hence we see that m(— 1) tends to zero for fi% — > but since it stays on the negative side, the matter era is unstable 
(one of the eigenvalues exhibits a positive divergence). So we can draw from Eq. (|74[) an important conclusion that 
the matter phase can only be obtained for fix = 0, i.e., the Starobinsky's (inflation) model previously discussed. 
In order to satisfy solar system constraints a particular version of this model was suggested with [30| 

H2 = 3 3/ Vi , and P = VSfil . (75) 

In that case, (f7^|) yields m(— 1) ~ 3.40 hence this case does not have a standard matter phase either. 
The model has two accelerated attractors: 

P 6 : (x 1 ,x 2 ,x 3 ) = (-2,3/2,3/2), Wcff = -2/3, , , 

Pi :(x 1 ,x 2 ,x 3 ) = (0,-1,2), w eff = -l. 

Thus depending upon the initial conditions the trajectories lie in the basin of attraction of either of these two points. 
This model corresponds to Class I. 
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Table I: Classification of f(R) dark energy models. 



8. m(r) = -0.2(1 + r)(3.2 + 0.8r + r 2 ) 

This model has been designed by hand to meet the condition for the Class IV. Note that this corresponds to 
the m(r) curve in the Class IV case shown in Fig. [H The corresponding f(R) Lagrangians are the solutions of the 
differential equation 

Rf.RR ( Rf.R 
: = m — 

f,R \ f 

which can be obtained numerically. This model obeys the conditions m(— 1) = and m' 5 > — 1 required for a saddle 
matter era P5 as well as the conditions (vo — l)/2 < me = 0.8 < 1 and m' 6 < —1 required for a stable accelerated 

point Pg in the region (D). The final accelerated attractor corresponds to the effective equation of state w e s ~ —0.935. 

p 1 

A model with similar properties but an analytical Lagrangian is f(R) = i?p rrT (i? + C) x -p (C ^ 0, p ^ 1) for which 
m(r) = — p(r + l) 2 /r, whose Pq intersection lies in the region (D) for 2 < p < 3.73 with w e s — ■ Here 

however the matter era has large eigenvalues so in fact is of very little duration and hardly realistic. A generalization 
to m(r) = —p(r + ro) 2 /r with ro slightly less than 1 works much better but then the Lagrangian is very complicated. 

9. Summary 

In Table I we summarize the classification of most f{R) dark energy models presented in this section. No model 
belongs to the Class IV except for the purposely designed cases given in the previous subsection, so we omit the Class 
IV column. The models which are classified in Class II at least satisfy the conditions to have a saddle matter era 
followed by a de-Sitter attractor. This includes models of the type / = R+aR~ n (— 1 < n < 0, a < 0), / = R(\ogaR) q 
(q > 0) and / = R + aR 2 — A. However this does not necessarily mean that these models are cosmologically viable, 
since it can happen that the matter era is too short or too long to be compatible with observations. In the next 
section we shall numerically study the cosmological viability of the above models. 

VI. SPECIFIC CASES: NUMERICAL RESULTS 

We will now use the equations derived in Section II in order to recover the cosmic history of given f(R) DE models 
and confirm and extend our analytical results. In all cases, we include radiation and give initial conditions at an 
epoch deep into the radiation epoch. As our aim is to check their cosmological viability, we tune the initial conditions 
in order to produce observationally acceptable values, namely 

On,o w 0.3 , r> rad , » 10- 4 . (78) 

In some cases we plot a 2-dimensional projection of the 3-dimensional phase space (xi,X2,xs) (no radiation) in 
Poincare coordinates, obtained by the transformation x' P ' = Xi/(l + d) where d = \J x\ + x\ + x^. 

A. f(R) = aR~ n 

Since m = — n — 1 in this case, the matter era is possible only when n is close to —1. So let us consider the 
cosmological evolution around n = — 1. As we already showed, the matter point P5 is stable for n < — 1. When 




f(R)=£KR°- 9 




Figure 4: Phase space in the plane (xi,X2) in Poincare coordinates for the model f(R) = aR 9 in the absence of radiation. 
Here and in the following plot, the dotted lines correspond to trajectories at the early stage, the continuous lines to those at 
the final stage. The circles represent critical points. The solutions approach either the <^MDE point P2 or the phantom point 
Pq. The point P5 is a saddle, but the trajectories do not approach this point if we take into account radiation. The point P3 
is an unstable node. 

n > —1, P5 is a saddle and both Pi and Pq are stable. In Fig. U] we show a 2-dimensional phase space plot for the 
model n = —0.9 in the absence of radiation. In fact the final attractors are either the <^MDE point Pi with w c s = 1/3 
or the phantom point Pq with w c g = —10.17. The point P5 with w e f[ = 1/9 is in fact a saddle point. However, if 
we start from realistic initial conditions around (x\, Xi, X3, 2:4) = (0,0,0,1) with the inclusion of radiation, we have 
numerically found that the trajectories directly approach final attractors (Pi or Pq) without reaching the vicinity 
of P5. Moreover as we choose the values of n closer to —1, the point P 5 becomes repulsive because of the positive 
divergence of an eigenvalue. These results show that the power-law models with n > —1 do not provide a prolonged 
matter era sandwiched by radiation and accelerated epochs in spite of the fact that the point P5 can be a saddle. 

B. f(R)=R + aR- n 

When n > one has m = —n — 1 < — 1 and m'(r) = n/r 2 > for P5 and Pq. In this case Pq is a stable attractor 
whereas P5 is a saddle. In the previous section we showed that the matter point P5 is disconnected to the accelerated 
point Pq since Pq exists in the region (A) . According to the results in Ref. we have only the following two cases: 
either (i) the matter era is replaced by the <^>MDE fixed point P2 which is followed by the accelerated attractor Pq, or 
(ii) a rapid transition from the radiation era to the accelerated attractor Pq without the </>MDE. Which trajectories 
are chosen depend upon the model parameters and initial conditions. In Fig. [5] we depict a 2-dimcnsional phase space 
plot for the model n = 1. This shows that the final attractor is in fact Pq and that whether the solutions temporally 
approach the saddle point Pi or not depends on initial conditions. 



f(R)=R+a/R 




Figure 5: Phase space projected on the plane (xi,X2) in Poincare coordinates for the model f(R) = R + a/R in the absence 
of radiation. For the initial conditions X2 > there are two solutions: either (i) the solutions directly approach the accelerated 
attractor P% or (ii) they first approach the saddle <^)MDE point P2 and then reach the attractor Pe- When X2 < initially, the 
trajectories moves toward X2 — > —00. Note that the point P3 is unstable. 
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Figure 6: The cosmic evolution of various quantities is shown for the model f(R) = R + a/R with a = — fi , fx/Ho = 11.04. 
The standard matter era is replaced by the 0MDE which corresponds to a oc t 1 ' 2 , w e g = 1/3 and fi m = 2. The redshift z a at 
which acceleration starts is z a = 0.4 and we have asymptotically in the future SIde = 1 and w c g — woe ~ —0.82 [see Eq. (|63|l] . 



In order to understand the evolution after the radiation era let us consider the model n = 1 without radiation. 
From Eqs. (jl]) and ([SJ we find that the evolution of the scale factor during the 0MDE is given by 

a(t) = (t/t i ) 1 / 2 + e(t)(t/t i ) 9 /\ (79) 
where the subscript V represents the value at the beginning of the <^>MDE. At first order in tit) we have 

« t ) = ,uH* m 1 (80) 



Notice that fi is of order Hq to realize the present acceleration. Since Hg <C Hi the parameter e(t) is in fact much 
smaller than unity. The scale factor evolves as a oc t 1 / 2 during the </>MDE, but this epoch ends when the second term 
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Figure 7: Projected phase space in Poincare coordinates for the model f(R) = R + aR 09 in the absence of radiation. The 
final attractor is the de-Sitter point Pi : (xi,X2, X3) = (0,-1,2). Note that both P2 and axe not stable unlike the model 
f(R) = aR 9 . In the left panel, a > 0: here P5 corresponds to m < with a large eigenvalue and therefore is unstable. In the 
right panel, a < 0: now the point P5 is a saddle with positive m, so it is possible to have a sequence of an oscillating matter 
phase followed by the late-time acceleration. We plot a single curve for clarity. 

in Eq. (|79|) gets larger than the zero-th order term. Hence the end of the 0MDE is characterized by 



After that the solutions approach the accelerated attractor P5. Equation (|81[) shows that the duration of the 0MDE 

depends on \i together with the initial conditions pm and Hi. The similar argument can be applied for any n < 
— l,7i > —3/4 with a correction growing as i 5 / 2 - 1 / 2 ("+ 1 ). i n Fig. El we plot the evolution of various quantities for 
n = 2. In this case the radiation era is followed by the </>MDE saddle point P2 with f2 m = 2 and w e s = 1/3. The final 
attractor is the accelerated point Pq with f^E = 1 and w e g = —0.82. As is clearly seen in the right panel of Fig. [S] 
we do not have a standard matter era with w e s = 0. 

Let us consider the case in which n is close to —1. When n < — 1 the point P5 is a stable spiral, so the matter era 
is not followed by an accelerated expansion as is similar to the power-law models. If n > —1, the de-Sitter point Pi 
is stable whereas the phantom point P§ is not. In Fig. [7] we show the phase space plot in a two-dimensional plane for 
n = —0.9. When a > 0, although the point P5 is a saddle, the solutions approach the attractor Pi without staying 
the region around the point P5 for a long time because m is negative. This tendency is more significant if n is chosen 
to be closer to —1, i.e. m — * —0. Hence one can not have a prolonged matter era in these cases as well. On the other 
hand, for a < 0, we have m — > +0 and there are oscillations around the matter era followed again by the attractor 
Pi. Then this latter case, belonging to Class II, can be cosmologically viable. 



When q > 0, we showed that the point P5 is a saddle for m(r < —1) > whereas both Pi and P§ are stable. In 
the previous section we showed that the only possibility is the trajectory from P5 to Pi. Hence the solutions starting 
from the radiation era reach the saddle matter point P5 first, which is followed by the de-Sitter point Pi. 

In order to obtain a prolonged matter period, the variables m and r need to be close to +0 and r = — 1, respectively, 
at the end of radiation era. If we integrate the autonomous equations with initial conditions r = X3/X2 — —1 (and 
smaller than —1) and X4 ~ 1, we find that the matter era is too long to be compatible with observations. In Fig.[S]we 
plot one example of such cosmological evolution for q = 1. This shows that a prolonged (quasi) matter era certainly 
exists prior to the late-time acceleration. The final attractor is the de-Sitter point Pi with w c tf = — 1. However in this 




(81) 



C. f(R) = R(logaR) q 



case the beginning of the matter epoch corresponds to the redshift z = l.lx 10 17 , which is much larger compared to 
the standard value z ~ 10 3 . The present value of the radiation energy fraction is fi ra d,o = 2.8 x 10~ 15 and is much 
smaller than the value given in Eq. (|78|) . 

This unusually long period of the matter era is associated with the fact that the point P5 is a saddle in the region 
r < — 1 but it is marginally stable in the limit r — > — 1 (i.e. m — * +0). Hence as we choose the initial values of r 
closer to —1, the duration of the matter period gets longer. In order to recover the present value of f2 ra d given in 
Eq. (|78|). we have to make the matter period shorter by appropriately choosing initial conditions at the end of the 
radiation era. In Fig.[TJ]we plot the cosmological evolution in the case where the end of the radiation era corresponds 
to z ~ 10 3 with present values r2 m o ~ 0.3 and f2 ra d,o ~ 10~ 4 . The energy fraction of the matter is not large enough 
to dominate the universe after the radiation epoch. Hence this case is not compatible with observations as well. 
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Figure 8: The cosmic evolution of various quantities for the model f(R) = RlogaR with initial conditions x\ — 10~ 5 , 
X2 = — 10~ 10 , £3 = 1.01 x 1CT 10 and X4 = 0.999 at the redshift z = 1.1 X 10 17 , corresponding to r = —1.01. In this case the 
matter era is too long relative to the standard cosmology. In fact the energy fraction of the radiation at the present epoch is 
n ra d,o = 2.8 x 10 -15 , which is much smaller than the standard value f2 ra d,o ~ 10 -4 . 
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Figure 9: The cosmic evolution of various quantities for the model f(R) = RlogaR with initial conditions xi = 10 -10 , 
X2 = — 10~ 7 , xg — 1.019 x 10 -7 and £4 = 0.999 at the redshift z = 3.15 x 10 6 , corresponding to r = —1.019. In this case we 
have f2 m ,o ~ 0.3 and f2 ra d,o ~ 10 -4 at the present epoch, but the matter era is practically absent. 



D. f(R) = Rexp(q/R) 

In this case the matter point P5 is a saddle, but one of the eigenvalues are 3 rather than close to 0. Numerically we 
find that the solutions do not reach the matter dominated epoch unlike the f(R) = R(\ogaR) q model with q > 0. In 
Fig.Hrjwe plot the cosmological evolution for this model corresponding to the present values f2 mj o ~ 0.3, ri ra d,o ~ 10 -4 . 
In this case the matter epoch is replaced by the ^MDE. It is possible to find a situation in which there exists a short 
period of the matter era, but we find that this case does not satisfy the conditions given by ([75)1 . Thus this model is 
not cosmologically viable in spite of the fact that it belongs to the Class II. 
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Figure 10: The cosmic evolution of various quantities for the model f(R) = Rexp(q/R) with initial conditions x\ — 0, 
X2 — 2.13 x 1CP 20 , xs — 5.33 x 10~ 21 and X4 = 0.99 at the redshift z — 3 x 10 5 . We see that the matter era is absent and is 
replaced by the 0MDE. 



E. m(r) = -0.2(1 + r)(3.2 + 0.8r + r 2 ) 

This model belongs to the Class IV, so the cosmological trajectories can be acceptable. In Fig.QT]we find that the 
matter epoch is in fact followed by a stable acceleration with w c s w —0.935. The transition between the various eras 
is not very sharp compared to the ACDM model, so it is of interest to investigate in more detail whether this model 
can be really compatible with observations. However this is beyond the scope of this paper. 




Figure 11: The cosmic evolution of various quantities for the model m(r) = —0.2(1 + r)(3.2 + 0.8r + r 2 ) with initial conditions 
xi = 10" 10 , x 2 = -10" 7 , x 3 = 1.000007 x 10~ 7 and x 4 = 0.999 at the redshift z = 3.5 x 10 6 . The model has an approximate 
matter dominated epoch followed by a non-phantom accelerated universe with u> e ff ~ —0.935. 



VII. CONCLUSIONS 

The f(R) dark energy models are interesting and quite popular attempts to explain the late-time acceleration. 
However it was recently found that the popular model f(R) = R + aR~ n with n > is unable to produce a matter 
era prior to the accelerated epoch [34| . In this paper we have attempted to clarify the conditions under which f(R) dark 
energy models are cosmologically viable. We first derived the autonomous equations pS jl - P^ which are applicable 
to general f(R) models. In Sec. Ill all fixed points are derived in such an autonomous system. By considering linear 
perturbations about the fixed points, we have studied their stabilities to understand the cosmological evolution in 
f(R) dark energy models. 

The main result of this paper is that we have identified four classes of f(R) models, depending on the existence 
of a standard or "wrong" matter era (</>MDE) and on the final acceleration. In practice, we have shown that the 
cosmology of f(R) models can be based on a study of the m(r) curves in the (r, m) plane and on its intersections 
with the critical line m = — r — 1. This provides an extremely simple method to investigate the cosmological viability 
of such models. In particular, we find that the Class I models correspond to the type of models in which the final 
acceleration is preceded by a so-called 0MDE phase characterized by a oc t 1 / 2 or in which the matter phase does not 



exist at all prior to the accelerated e poc h. These models are clearly ruled out, e.g. by the angular diameter distance 
of the CMB acoustic peaks, see Ref. [34]. This is by far the largest class and only a few special cases belong to other 
three. 

The general conditions for a successful f(R) model can be summarized as follows: 

• A f(R) model has a standard matter dominated epoch only if it satisfies the conditions 

m(r) w +0 and m'(r) > -1 at r ss -1 , (82) 
where the second condition is required to leave the matter era for the late-time acceleration. 

• The matter epoch is followed by a de-Sitter acceleration (w e g = — 1) only if 

< m(r) < 1 at r = -2 or m(r) = -r - 1 ->■ ±oo (Class II) . (83) 

• The matter epoch is followed by a non-phantom accelerated attractor (w c g > — 1) only if m = —r — 1 and 

(V3 - l)/2 < m(r) < 1 and m'(r) < -1 (Class IV) . (84) 

Moreover, the curve m(r) must connect with continuity the vicinity of the matter point P5 : (r, m) — (—1, 0) with one 
of the accelerated regions. The Class II and IV models are characterized by m(r) curves that satisfy these requirements 
and lead therefore to an acceptable cosmology. 

In the Class III models the curve m(r) intersects the critical line at m small and negative. In this case the saddle 
eigenvalue takes a very large real value and the matter era is practically unstable and therefore generically very short. 
Moreover, most trajectories will be attracted by the strongly phantom attractor with w c g < —7.6 which is in contrast 
with observations. 

The cases with m'(r) = — 1 or m'(r) — > ±00 at the critical points are not covered in our linear approach and a higher- 
order or numerical analysis is necessary. Also the power-law model f(R) = R~ n is a rather special case in a sense that 
it gives a transition from the quasi matter era to the strongly phantom epoch with a constant negative m. However we 
showed that this model is not cosmologically acceptable because of the absence of the prolonged matter epoch. We have 
also studied analytically and numerically models like f(R) = R + atR~ n , R p (\og aR) q , R p exp(qR), RP exp(q/R), (R a — 
A) b and others and have confirmed the conclusions drawn from the m(r) approach. See Table I for the summary of 
the classification of a sample of f(R) dark energy models. 

As we have seen, the variable m = Rf,Rn/f,R plays a central role to determine the cosmological viability of f(R) 
dark energy models. The ACDM model, f(R) = R — A, corresponds to m = at all times, which thus satisfies the 
condition for the existence of the matter era (m ss 0) followed by the de-Sitter point at m(r = —2) = 0. The difference 
from the line m = characterises the deviation from the ACDM model. If the devaition from m = is small, it is 
expected that such models are cosmologically viable. 

We conclude with a comment concerning a possible signature of f(R) cosmology. The standard matter era can 
be realized with m — * ±0. As we have seen, in all successful cases we analyzed in this work, the matter era is 
realized through damped oscillations with positive m. This raises the obvious question of whether such oscillations 
are observable and whether they could be taken as a signature of modified gravity. This question is left to future work. 
An additional interesting direction to investigate is the evolution of cosmological perturbations in f(R) dark energy 
models in order to confront with the datasets of CMB and large scale structure along the lines of Refs. [H, |45|, [4(| . 
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